Importance of electronic self-consistency in the TDDFT based treatment of 

nonadiabatic molecular dynamics 
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A mixed quantum-classical approach to simulate the coupled dynamics of electrons and nuclei 
in nanoscale molecular systems is presented. The method relies on a second order expansion of 
the Lagrangian in time-dependent density functional theory (TDDFT) around a suitable reference 
density. We show that the inclusion of the second order term renders the method a self-consistent 
scheme and improves the calculated optical spectra of molecules by a proper treatment of the coupled 
response. In the application to ion-fuUerene collisions, the inclusion of self-consistency is found to be 
crucial for a correct description of the charge transfer between projectile and target. For a model of 
the photoreceptor in retinal proteins, nonadiabatic molecular dynamics simulations are performed 
and reveal problems of TDDFT in the prediction of intra-molecular charge transfer excitations. 



I. INTRODUCTION 

Beginning with the work of Zangwill and Soven 
the generalization of density functional theory to time 
dependent phenomena (TDDFT) has become an impor- 
tant tool in the description of laser-matter interaction. 
The possible applications are diverse and range from the 
calculation of spectra (linear optical 0, 0; Q|j circular 
dichroism ^sl, 'gI, resonant Raman '^) to the evaluation 
of properties (polarization 8, 9], hyperpolarization UOl) 
up to studies of high harmonic generation [ill ll2L ll.^] 
and photochemical reactions jljj . The formal justifica- 
tion of TDDFT was laid by Runge and Gross [13 , who 
showed that the exact many body electron density can be 
obtained from single-particle mean field equations. The 
solution of these time dependent Kohn-Sham (TDKS) 
equations can be obtained either pcrturbatively in the 
small amplitude limit 16| or by direct numerical integra- 
tion in the time domain |jl7i ] . Both approaches have their 
inherent merits and disadvantages. 

In the linear response regime, for example, the problem 
can be recast in an eigenvalue equation in the particle- 
hole representation. This allows for an interpretation of 
optical spectra in terms of contributing single-particle 
transitions and also for a symmetry assignment of the 
states. Moreover, transitions with vanishing oscillator 
strength can be located, like dark singlet or generally 
triplet states [sj . One of the drawbacks of this approach 
is the rather poor numerical performance with a scaling 
of , where N is the number of electrons. It should be 
noted, however, that the CPU time as well as the mem- 
ory demand can be significantly reduced when iterative 
procedures like the Davidson algorithm are employed. 

In terms of scaling behavior, the numerical integration 
of the TDKS equations is much more favorable. Here, 
only the set of occupied orbitals needs to be treated. Be- 
cause the propagation involves only matrix-vector prod- 
ucts, linear scaling can be achieved for large systems [T^ . 
Moreover, since this approach is not restricted to small 
intensities, non-linear effects like harmonic generation or 



multiphoton processes can be addressed. Another advan- 
tage of working in the real time domain is the possibility 
to study systematically the effect of different pulse shapes 
of the laser field on observables such as ionization |20l| . 
With todays femtosecond laser sources, this is currently 
an active field of experimental research |2ll |. 

Nearly all first principles applications of the real time 
approach have been limited to systems with fixed nu- 
clei. Clearly, it would be highly desirable to study the 
motion of the coupled system of electrons and nuclei, 
which would allow one to address problems like laser in- 
duced vibrational excitation or photochemical reactions. 
Since the time step of such molecular dynamics simu- 
lations is set to attoseconds by the ultrafast electronic 
motion, only small systems with a few degrees of free- 
dom can be treated in an ab initio frame work 0, ■ 
Consequently, approximate TDDFT methods are quite 
successful in this domain of application. 

Different groups contributed to this field and used their 
developments in a variety of different studies 13, 
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 3g. 
In all these approximate schemes only the valence elec- 
trons are treated explicitely and the TDKS orbitals are 
expanded in a limited (usually minimal) basis of atomic 
orbitals. The Lagrangian, which is a functional of the 
time dependent density, is then expanded around a static 
reference density up to a certain order. In zeroth order 
the Hamiltonian depends only on the reference density, 
which permits the calculation of the necessary matrix el- 
ements once and for all. In this respect, the methods are 
similar to tight-binding approaches, although no fitting 
to experimental data is performed. 

The purpose of this work is to analyze the implica- 
tions of extending the mentioned expansion, since all 
studies so far were restricted to zeroth order. After a 
more detailed description of the problem in Sec. ^ we 
test the extension in the determination of optical spec- 
tra in Sec. IIII Al as well as for nonadiabatic molecular 
dynamics in Sec. IIII 51 Finally, we perform an investiga- 
tion of the photochemical reaction of a retinal analogue. 
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a chromophore which exhibits an ultrafast radiationless 
deactivation in nature. These apphcations in quite differ- 
ent areas of molecular physics are intended to investigate 
the transferability of the method and also to illustrate 
the possibilities offered by an approximate solution of 
the TDDFT equations. 



II. METHOD 

In order to study the dynamics of a coupled system 
of electrons and nuclei, the equations of motion (EOM) 
need to be determined. While the electronic EOM in 
the framework of DFT is given by the well known time 
dependent Kohn-Sham equations, the nuclear EOM or 
force equation is not a priori evident. It can be derived 
either by exploiting the fact that the total energy is a 
conserved quantity, or by applying the Lagrange formal- 
ism. We follow the latter approach here and define the 
following Lagrangian, which depends on the TDKS states 
'$i{rt) and the nuclear positions R^i: 



(1) 



(M/.(r<)| H[p]{rt) - I- |*,(rt)) - E^c - E, 



with p{Yt) — Here the first term is the 

classical kinetic energy of the ions, while the remaining 
terms in Eq. can be obtained from the TDDFT ac- 
tion functional under the assumption that the exchange- 
correlation (xc) contributions are local in time |44j . In 
this widely used adiabatic local density approximation, 
standard ground state functionals can also be used in the 
time dependent context simply by evaluation at the time 
dependent density. Thus, the Hamiltonian H[p]{vt) in 
Eq. takes the common DFT form. Furthermore, ii^DC 
represents the double counting terms 



' p{vt)p{v't) 



Exc[p] - / Vxc[p]p{vt), 



and Eii the ion-ion repulsion (Here and in the following 
J dr' is abbreviated as /', and J dr as /). 

We now proceed by applying the same kind of approx- 
imations as were used in the derivation of the density 
functional theory based tight-binding (DFTB) method 
[4(ll l4ll| from static DFT. To k eep the presentation con- 
cise we refer to some reviews |38l l39l |. which provide a 
more detailed description of the basic concepts, practical 
realization and accuracy of the ground state DFTB ap- 
proach. Here, we only report aspects, which are specific 
for the generalization to the time-dependent case. In a 
first step, the Lagrangian is expanded around a reference 
density po(r), p(rO = Po{^) + Sp{rt), which is given as a 
superposition of atomic ( grou nd state) densities. In con- 
trast to our earlier work |39j, we now include terms up 



to second order in the density fluctuations Sp{rt): 
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Please note that in this expansion all contributions which 
are linear in Sp are captured by the second term in 
Eq. ^ through the TDKS states. The terms in Eq. (|3b|) 
can now be subsumed as -Erep i a- sum of short ranged pair 
potentials, which depend only on the atomic species and 
the interatomic distance. Since i?icp is a functional of the 
time independent reference density po only, it is exactly 
the same as used in the ground state DFTB scheme. The 
second order term of Eq. l|llc|) , which is the focus of this 
work, is approximated as follows: 
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Here the AqA{t) denote atomic net MuUiken charges 
AqAit) = qA{t)~qT''^'°"' (5) 
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where the coefficients h^i(t) are defined by the expansion 
of the TDKS states in a basis of non-orthogonal atomic 
orbitals <j)^j,{v — R^) 



(6) 



Further, 



which build the overlap matrix S*^^ = (0/^1 ( 
the function ^ab in Eq. 0] interpolates between a pure 
Coulomb interaction for large interatomic distances and 
an element-specific constant in the atomic limit. This 
numerically evaluated parameter includes the effects of 
exchange and correlation and is directly related to the 
chemical hardness of the atomic species Taking 
the coefficients and nuclear positions as generalized coor- 
dinates, the evaluation of the Euler-Lagrange equations 
leads to the desired equations of motion. The electronic 
motion obeys: 
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In zeroth order the Hamiltonian reduces to the first term 
in Eq. I|7bp and depends solely on the reference den- 
sity po- For systems in the excited state or in charged 
or hetcronuclcar structures, the electronic density dif- 
fers significantly from a simple superposition of atomic 
ground state densities. To a certain extent this differ- 
ence is already captured at the zeroth order level, since 
the coefficients which solve Eq. ()7a|l correspond to a time- 
dependent density different from po- This is similar to 
the situation in empirical tight-binding schemes for the 
ground state, where even certain ionic crystals are suff- 
ciently well described j43|. However, a consideration of 
the full Hamiltonian in Eq. I)7b|l leads obviously to a more 
balanced treatment, because dynamical changes in the 
electron density are explicitly included in a selfconsis- 
tent fashion. In this way, one can even hope to correctly 
describe the large amplitude motion induced by intense 
laser fields. 

Solution of Eq. (|7a|l requires an iterative procedure 
with timesteps in the attosecond regime. A symplectic 
algorithm is used for this task [s^l , which is based on the 
Cayley representation of the time evolution operator and 
conserves the norm of the wavefunction exactly. Besides 
the Hamiltonian and overlap matrices (see Ref. [s^ for 
details of the construction) , Eq. (|7a|l contains also the 
nonadiabatic coupling matrix {<Pii\g^4'u} , in which all 
on-site elements {(j}fj,,4>v on the same atom) are set to zero 
. This allows one to relate the remaining elements to 
a simple derivative of the overlap matrix. 

Variation of the Lagrangian with respect to the nu- 
clear coordinates leads to the following expression for the 
forces: 
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FIG. 1: Schematic illustration of trans-butadiene {dHe 



Todorov pointed out, that for an incomplete basis the 
force equation has to be augmented by additional velocity 
dependent terms, which should become important in high 
energy collisions. Interestingly, omission of these terms 
does not violate energy but only momentum conserva- 
tion. Hence, it is useful to monitor the total momentum 
(electronic -I- ionic) of the system if Eq. © is used, as in 
every practical calculation the basis set is incomplete. 

Finally, in order to simulate the interaction with elec- 
tromagnetic fields, the vector potential A(rt) needs to 
be incorporated, which is done via minimal coupling, 
p — |A. A numerically efficient approximation was pro- 
posed by Graf and Vogl 45] and later by Allen and 
is given by: 



i?^^[r, p, A(t)] = exp 



he 



{Ra - Rb)A(0 



(8) 



xi/,,^[r,p]; pe^,i/GB, (9) 

which relates the desired field dependent Hamiltonian to 
the already known matrix elements of the unperturbed 
one. Expression @ was derived under the assumption 
that the radiation wavelength is much larger than the 
molecular system under study, which is usually well ful- 
filled for frequencies in the optical range. It should be 
mentioned, that Eq. @ in principle holds for arbitrarily 
strong fields in contrast to the electric dipole approxima- 
tion. 

If the interest is just on calculating optical spectra, 
rather than the molecular motion initiated by a laser 
pulse, only Eq. H7a() needs to be solved for a fixed ge- 
ometry. Following the approach of Yabana and Bertsch 
[r7. 46] the field in Eq. ^ is turned on only at a certain 
instant of time, which populates the complete manifold of 
excited states. The time dependent dipole moment d{t) 
can then be used to calculate the dynamic polarizability: 



Since the Hamiltonian is time dependent due to an ex- 
ternal field or nuclear motion, the molecular orbital coef- 
ficients will in general represent a coherent superposition 
of different eigenstates of the system. In this case, the nu- 
clei move in a mean potential according to Eq. |S1 rather 
than being restricted to a particular Born-Oppenheimer 
surface as in conventional adiabatic MD approaches. In 
fact, due to the coupling of the EOM, energy can be freely 
exchanged between the electronic and ionic subsystems 
as long as the total energy of the system is conserved. 
Equations of motion that are equivalent to the ones re- 
ported h ere, have been derived earlier by Saalmann and 
Schmidt (23 as well as Todorov 0. Including the sec- 
ond order correction, they are solved here for the first 
time in an actual calculation. 



{d{t) - d(0)) dt, 



(10) 



and the dipole strength function S{uo) = 2uo /ir'^aluo) of 
the system; a quantity which can be directly compared 
to experimental spectra. 



III. APPLICATIONS 

A. Optical spectrum of trans-butadiene 

As a first application of our method, we examine a pro- 
totypical TT-system, trans-butadiene (Fig.0. The optical 
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proach, excited state singlet energies uji are given by the 
solution of the following eigenvalue problem: 
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FIG. 2: Dipole strength of trans-butadiene as given by the 
time-dependent DFTB method in zeroth and second order 
• Shown is an average over different molecular orienta- 
tions with respect to the polarization of the vector potential. 
Results obtained in the DFTB linear-response implementa- 
tion of TDDFT are shown as stick spectrum. The associated 
oscillator strengths have been normalized to the maximum of 
the highest energy peak. 



spectrum of this molecule has been the subject of numer- 
ous quantum chemical studies (see e.g. Ref. H and ref- 
erences therein). Recently, also a detailed investigation 
of the molecular dynamics in the excited state appeared 
[T^ . Dou et al. employed the DFTB method described 
in this work without the second order correction. Our 
interest here is to analyze the implications of including 
this term. To this end, we first relaxed the molecule with 
the ground state DFTB method and recorded the optical 
spectrum according to the prescription given in Sec. ITU 
After applying a vector potential of A = 0.0125 gauss cm, 
the Kohn-Sham orbitals were propagated for 38.7 fs with 
a time step of 12 as. Since the finite sampling introduces 
spurious negative parts in the imaginary part of the po- 
larizability, the dipole moment was damped with a factor 
of e-*^* (k = 0.3 eY/h), hkc in Ref. |l3 This also sim- 
ulates dephasing or other line broadening effects which 
would appear more naturally in a more complete theory. 
The resulting spectrum with and without second order 
correction is shown in Fig. |21 As can be seen, the max- 
imum absorption in the latter method is located at 4.21 
eV. This is exactly the difference of the LUMO (lowest 
unoccupied molecular orbital) and the HOMO (highest 
occupied molecular orbital) energies of the ground state 
DFTB method. If the second order term is included in 
the calculation, the absorption is strongly blue shifted to 
5.56 eV, which is in good agreement with the experimen- 
tal value of 5.8 eV 0|. Along with the energy shift, a 
reduction of absorption strength is also observed. 

To better understand the origin of these changes, we 
also performed calculations with our implementation of 
the TDDFT linear response formalism 50]. In this ap- 
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[uj^^S^kSji + 2^/uJ~KijM^/^i\ Hi = ^1 Fij- (11) 



Here the Uij are energy differences between unoccupied 
orbitals j and occupied orbitals i, while the so called cou- 
pling matrix K describes the change of the SCF potential 
due to the induced density. As Eq. shows, the effect 
of the coupling matrix is not only to shift the true excited 
state away from simple orbital energy differences, but 
also to couple different single-particle transitions. The 
explicit form of K is given by: 



ij.kl 



?/'i(r)V'j(r) 

Spiv') 



(12) 



r — r' 



which is nothing else than the second order term of 
Eq. H3c|l . when the induced density is expanded in 
particle-hole states. The results of the linear response 
calculations with and without the coupling matrix con- 
tribution are given as stick spectrum in Fig. ^ Obvi- 
ously, there is a perfect match beween the real time and 
linear response approaches to TDDFT both in the ener- 
getical position of the states and the oscillator strength. 
This equivalence had to be expected, since the linear 
response approach amounts to a perturbative solution 
of the TDDFT equations in the small amplitude limit. 
However, to our knowledge this has so far never been 
shown in practical calculations. 



B. — Ceo collisions 

Ion-cluster collisions provide an ideal application field 
for approximate TDDFT molecular dynamics simula- 
tions. This is because the number of degrees of free- 
dom is usually too large to be treated with first princi- 
ples calculations and also nonadiabatic effects are strong 
and important. Depending on the velocity of the pro- 
jectile, collisions can induce vibrational or a combina- 
tion of electronic and vibrational excitations of the clus- 
ter, where the latter type cannot be described in con- 
ventional Born-Oppenheimer dynamics. In this context, 
Kunert and Schmidt undertook a systematic investiga- 
tion of ion-fuUerene collisions and provided an explana- 
tion for seemingly conflicting experimental observations 
|34j . Their nonadiabatic quantum molecular dynamics 
(NA-QMD) method is essentially equivalent to the DFTB 
scheme described here in the zeroth order approximation. 
Accordingly, it is interesting to see whether the second 
order correction is of any benefit in these kind of simula- 
tions. 

For the special case of 0+ — Cgo collisions we performed 
calculations for different values of the impact velocity {v 
= 0.01 . . . 0.45 a.u.) and impact parameter {b =2.0,7.5 
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FIG. 3: Total kinetic energy loss AE of the projectile in the 
center of mass system for an impact parameter of b = 2.0 a.u. 
For each velocity three trajectories with random cage orien- 
tation were calculated; error bars correspond to the standard 
deviation. Dark squares: second order, open circles: zeroth 
order DFTB results, open triangles: loss due to electronic ex- 
citation, obtained from the difference of time-dependent and 
ground state energy in the zeroth order DFTB scheme. 



a.u.) for randomly oriented fuUerene cages. Special care 
had to be taken in the definition of the initial conditions 
of the EOM, since the C"*" — Ceo configuration does not 
correspond to the ground state of the system. For that 
reason, we performed separate ground state calculations 
for the two subsystems and combined the resulting KS 
orbitals to obtain the desired charge state. The initial 
ion-cluster distance was chosen large enough to prevent 
any interaction and the system was then left to evolve 
freely according to the EOM [Eq. |(7a|) and ©]. 

Fig. O depicts the total kinetic energy loss /S.E in the 
center of mass system that the projectile experiences due 
to the collision. It can be directly compared to the re- 
sults in Fig. 3 of Ref. |3 which were obtained for a 
single fixed collision geometry. As already shown there, 
the vibrational excitation of the fullerene dominates for 
smaller velocities {v < 0.1 a.u.), while mostly electronic 
excitation is responsible for the energy loss at larger im- 
pact energies. As the impact velocity increases, AE first 
rises, peaks around 0.05 a.u and shows a weak increase 
beyond 0.1 a.u. in our calculations. This is in variance 
with the results of Kunert and Schmidt [sSl , which claim 
velocity-independent excitation energies in the high ve- 
locity range. Taking the strong dependence on the col- 
lision geometry into account, an extensive phase space 
sampling would be necessary to resolve this issue, which 
is outside the scope of this work. 

Turning now to a comparison of the predictions of 
DFTB with and without second order correction, we find 
only marginal differences in the results of both methods. 
Such a difference could have been expected for high im- 
pact velocities, but with such a large amount of energy 



FIG. 4; Charge of the carbon atom after the collision with 
Ceo for different values of the impact velocity and an impact 
parameter of b = 7.5 a.u. . For each velocity three trajecto- 
ries with random cage orientation were calculated; error bars 
correspond to the standard deviation. Dark squares: second 
order, open circles: zeroth order DFTB results. 



deposited in the cluster, finer details of the electronic 
structure seem to have negligible influence. 

The second order correction has however implications 
for other observables. Fig.^shows the charge of the pro- 
jectile after the collision. Here fractional charges need 
to be understood in the probabilistic interpretation of 
quantum mechanics, since in the simulations the sys- 
tem remains in a superposition of eigenstates with in- 
teger charges also asymptotically. For higher velocities, 
which correspond to higher electronic excitation as men- 
tioned above, the charge given by the zeroth order DFTB 
method is significantly more negative than the second or- 
der one. This can be explained by a larger contribution 
of the asymptotic — Cqq state, which in the zeroth 
order approximation is located only slightly higher in en- 
ergy as the initial C"*" — Cgg state. As the results in Tab.|ll 
show, this energy ordering is in striking contrast with the 
one given by the second order DFTB method and ex- 
periment. Although calculations of charge transfer cross 
sections have been performed in a zeroth order scheme 
'35,'335, the results of this section suggest, that in general 
a more advanced treatment is absolutely necessary. 



C. Protonated Schiff base photodynamics 

As a last application we study the photodynamics of 
the retinal molecule, which is of special importance in the 
field of biology. This chromophore is found in a variety 
of proteins, where it initiates quite different reactions in 
the cell. In bacteriorhodopsin and halorhodopsin for ex- 
ample, light absorption of retinal triggers the membrane 
transport of protons and chloride ions, respectively. In 
contrast, it starts a cascade of reactions that initiate the 
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TABLE I: Energy ordering of different asymptotic states of 
the singly positively charged C — Ceo system in eV. The DFTB 
energies with and without second order correction were ob- 
tained by separate calculations of the two subsystems and ad- 
dition of the results. For Ceo, all calculations were performed 
at the DFTB optimized geometry of the neutral species. The 
experimental results were obtained from measured ionization 
potentials and electron affinities 51, 52, 53]. The most stable 
state of each method was set to zero energy. 
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FIG. 5: Top: Structure of 11-cis retinal, which is found in 
dark adapted rhodopsin and transforms to the all-trans form 
upon absorption of light. Bottom: Protonated Schiff base 
model used in this study. 

vision process in rhodopsin. In all these systems, the 
retinal is known to isomerize around a specific double 
bond. The quantum yield of the photoreaction is partic- 
ularly high ($ w 0.7) and the deexcitation to t he g round 
state occurs in no more than 200-500 fs 1551 . Be- 
cause of these unusual features, this system provides an 
interesting subject for a theoretical investigation. For 
a complete understanding of the retinal photodynamics 
it would certainly be necessary to include the full pro- 
tein environment in such an investigation. However, im- 
portant information can already be drawn from the ex- 
amination of small retinal analogues like the protonated 
Schiff bases (Fig.O. These models share a polyene chain 
of alternating single and double bonds with retinal, as 
well as the positively charged NH^ Schiff base group, 
which is crucial for the function of the chromophore in 
the protein. High level quantum chemical CASPT2 cal- 
culations on these analogues revealed, that after absorp- 
tion the system moves out of the Franck-Condon region 
along the C=C stretch normal coordinate (see Ref. 
and references therein) . After inversion of single and dou- 
ble bonds, torsion around the central C=C bond sets in. 
A barrierless path then leads to a conical intersection of 
ground and excited state, where efficient deactivation to 
the photoproduct occurs. In line with Stark spectroscopy 



[57I Is^ l , C ASPT2 theory predicts a large charge transfer 
from the Schiff base group to the other terminus upon 
excitation, that increases along the excited state path- 
way. 

Recently, we investigated the excited state potential 
energy surface (PES) obtained from static DFTB calcu- 
lations in the linear response formalism and found se- 
vere deviations from the CASPT2 results |^. In fact, 
the only barrierless paths found to conical intersections 
with the ground state involved single rather than double 
bond isomerization (in accordance with ab initio TDDFT 
calculations). An intrinsic reaction coordinate connect- 
ing Franck-Condon point and correct conical intersection 
(as described in the CASSCF model) includes a signif- 
icant barrier. Thus an efficient and ultrafast reaction 
seems to be unlikely at the DFT level of theory. How- 
ever, one should keep in mind, that at finite temperature 
molecules posses a significant amount of kinetic energy 
already at the Franck-Condon point, which allows the 
system to sample a large fraction of the excited state 
PES. Hence, the minimum energy path might not neces- 
sarily provide a representative description of the photo- 
dynamical pathway. Moreover, nonadiabatic transitions 
are not restricted to conical intersections. They can also 
occur in regions where there is a finite gap between the 
ground and excited state surfaces, especially when the 
atomic velocity is high. Our interest is therefore to per- 
form nonadiabatic molecular dynamics simulations of the 
PSB model system to see whether the discrepancies be- 
tween DFT predictions and experiment remain at a full 
dynamical level. 

A complete description of the photochemical process 
would in principle require a full phase space sampling 
prior to excitation. Since the maximum absorption is 
strongly geometry dependent, we nevertheless take only 
the relaxed geometry of the ground state minimum into 
account. At the Franck-Condon point random velocities 
corresponding to a temperature of 300 K are assigned to 
the atoms and the system is left to evolve freely without 
further constraints. The excitation itself is induced by 
a Gaussian shaped laser pulse with a central frequency 
of 3.45 eV, that is slightly detuned with respect to the 
maximum absorption to avoid population of higher ex- 
cited states close in energy. The maximum absorption 
itself is located at 3.88 eV in the second order DFTB 
method and agrees well with first principles TDDFT cal- 
culations m (4.03 eV) as well as CASPT2 ^] (4.02 eV) 
and experiment (3.85 - 4.25 eV HJ). Similar to the case 
of trans-butadiene, the excitation energy is strongly un- 
derestimated at 2.66 eV, if the second order correction 
is not taken into account. Since the nonadiabatic excita- 
tion process depends strongly on the gap between ground 
and excited state PES, the DFTB method without elec- 
tronic selfconsistency will therefore provide an unrealistic 
description of the photodynamics and will not be used 
here. The duration and fiuence of the laser pulse were 
chosen to be 5.9 fs and 5.4 mJ/cm^ (A = 0.7 gauss cm). 
For these parameters, the total energy of the system rises 
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FIG. 6: Bond alternation in the PSB model system for all 
calculated trajectories, estimated as the mean bond length 
difference between neighboring C-C single and double bonds. 



50 
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FIG. 7: Torsion angle around the central double bond of the 
PSB model for all calculated trajectories. 



by 3.88 eV, i.e. a one-photon transition to the first ex- 
cited state is simulated. For a total simulation time of 
1.1 ps, 100 trajectories were propagated with a timestep 
of 12.1 as. This guaranteed an energy conservation of 
AE/E w 10~^. With these parameters, one trajectory 
took about 22 minutes CPU time on an Intel Xeon 3.06 
GHz processor. 

The results of the simulations show that the initial dy- 
namics on the excited PES are dominated by C-C stretch 
motions with large amplitudes up to 0.15 A. This is in 
agreement with resonant Raman studies on the PSB in 
solution and in the rhodopsin protein Ib^ . l65l |. In 
contrast to the mentioned CASPT2 calculations however, 
the bond alternation is only inverted for roughly half of 
the trajectories as shown in Fig. Moreover, the time- 
averaged bond alternation of 0.060 Ais only slightly re- 
duced with respect to the ground state minimum (0.063 
A) . Considering now the dihedral angle which represents 
the torsion around the central double bond. Fig. [3 shows 
that none of the 100 trajectories resulted in a success- 
ful isomerization within 1 ps. We find much larger am- 
plitudes for the rotation around single bonds, with tor- 
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FIG. 8: Excitation energy of the PSB model versus time in 
eV. The inset shows the transition to the excited state due to 
the applied laser pulse. 



sion angles up to 80° for certain trajectories, although 
also here no isomerization is completed in the simulation 
time. This preference of single over double bond iso- 
merization in DFT based methods was already found in 
the static investigation of Ref. Another prediction of 
CASPT2 theory, which is in nice agreement with Stark 
spectroscopy has already been mentioned. The Sq — Si 
excitation and the subsequent motion on the Si PES is 
known to involve a large charge transfer away from the 
Schiff base group. We however, do not observe any sig- 
nificant charge transfer in our simulations. 

Finally, it is interesting to analyze whether decxcita- 
tion to the ground state occurred. In Fig.|Slthe excitation 
energy of the system is shown, which is given by the dif- 
ference of the time-dependent and ground state energy 
at the same geometry. Directly after the end of the laser 
pulse, a small reduction of the excitation energy is ob- 
served (~ 0.3 eV), which is related to an elongation of all 
bonds in the PSB model. There is little change from this 
point on and nonadiabatic transitions, which would man- 
ifest in a sharp drop of the excitation energy, do not oc- 
cur. In our DFT based treatment, deactivation is there- 
fore predicted to happen on longer timescales, presum- 
ably involving the slower processes of internal conversion 
and fluorescence. At any rate, ultrafast isomerization 
with a concomitant intersection of ground and excited 
state is not found to be the dominant process. This is 
in stark contrast to the results of Vreven et al. , who 
showed that double bond isomerization can occur in less 
than 100 fs for the model at hand. 

The simulations of this section could be extended by 
computing a larger number of trajectories or a longer 
propagation time. Considering the narrow distribution 
of the results presented in Fig. |H1 to |H1 it is not very 
likely that an improved sampling would reveal new in- 
formation. Extension of the simulation time might seem 
advantageous in light of the experiments by Logunov et 
al. |6^, who measured an excited state life time of 2-3 
ps for the full retinal chromophore [Fig. [S] (top)] in solu- 
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tion. However, for the shorter PSB model used in this 
study, the excited state PES has been found to be much 
steeper |5fill in hne with the mentioned study of Vreven 
et al. [63 ■ Hence, the photochemical process should be 
completed in the chosen simulation time of 1 ps. 

To summarize, we find that our simulations disagree 
in most aspects with CASPT2 results and experiment. 
At the same time we confirmed the static investiga- 
tions of the DFT based potential energy surfaces from 
Ref . Is^ bv dynamical calculations. Given that, (i) time- 
dependent DFTB and first principles TDDFT, as weU 
as, (ii) TDDFT with different exchange-correlation func- 
tional (local, gradient corrected or hybrid) yield qualita- 
tively the same picture po] , a correct DFT based descrip- 
tion of the retinal photodynamics is highly unlikely. Only 
very recently, failures of TDDFT in the description of 
inter- molecular charge transfer states were reported |68j . 
Exchange-correlation functional that address this short- 
coming have also been recently proposed JQ,, 71]. It 
will be interesting to see whether these developments can 
also remedy the problems with intra-molecular charge 
transfer found here. 



IV. SUMMARY 

In this work, we presented a mixed quantum-classical 
approach to simulate the coupled dynamics of electrons 
and nuclei. The method is based on a second order ex- 
pansion of the TDDFT Lagrangian around a suitable 
reference density. We showed that the inclusion of the 
second order term improves both qualitatively and quan- 
titatively the optical spectrum of molecules. For trans- 
butadiene a strong blueshift of the absorption was ob- 
served together with a significant reduction of oscillator 
strength. In this context, the analogy with the linear 
response approach to TDDFT revealed that the zeroth 
order treatment of the Lagrangian corresponds to an un- 



coupled response that neglects collective effects. More- 
over, experience with the linear response implementation 
suggests that this absorption shift is quite general and 
especially large for n — tt* transitions. For n — ir* excita- 
tions however, the coupling is usually weak and realistic 
results might already be achieved in the zeroth order ap- 
proximation. In the simulations of high energy collision 
of Sec. IIIIBI there is little difference in the predictions 
of both schemes when the interest is in energy transfer 
only. This is because the process is dominated by vibra- 
tional rather than electronic excitation in the regime of 
low energy transfer, where the different level structure 
could be resolved. Considering now the charge transfer, 
we found important differences for the C+ — Ceo system, 
which were attributed to the problematic description of 
the asymptotic states in the zeroth order scheme. Charge 
transfer played also an important role in the nonadiabatic 
molecular dynamics simulations of the protonated Schiff 
base. In contrast to experiment we found no isomeriza- 
tion. This result is a negative, but we think important 
one. It should be stressed, that this failure is not in- 
troduced by our approximations but already inherent in 
TDDFT itself. 
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